##########################################################################
# HANGARTNER, LAUDERDALE, SPIRIG (2025) INFERRING INDIVIDUAL PREFERENCES #
##########################################################################

################################################
########### RESULTS, TABLES, FIGURES ###########
################################################


rm(list = ls())

stringsAsFactors <- FALSE
options(scipen = 999)

#############
# Libraries #
#############

library(dplyr)
library(rstan)
library(stargazer)
library(xtable)

#############
# Functions #
#############
 
round3 <- function(x) formatC(x,digits=3,format="f")
find_judge <- function(x) sum(grepl(a, x))

########
# DATA #
########

judgedataloc <- "Data/jdata.Rda"
load(judgedataloc)

posteriorfilelocation <- "SavedPosterior.Rdata"
load(posteriorfilelocation)

materieloc <- "Data/materie_codes.Rda"
load(materieloc)


##############################################################
## Do some judges systematically get weaker/stronger cases? ##
##############################################################

Data_test <- Data[!is.na(Data$lawyer),]
nrow(Data_test)

# add materie text
Data_test$MaterieAggregated <- materie_agg_codes$Text_en[match(Data_test$MaterieAggregated,as.numeric( materie_agg_codes$Code))]
Data_test$MaterieAggregated <- gsub("Reconsideration Request", "RR", Data_test$MaterieAggregated)
Data_test$MaterieAggregated <- gsub("Inadmissibility of Asylum Request", "Inadmissibility of Request", Data_test$MaterieAggregated)
Data_test$MaterieAggregated <- gsub("Return and Enforcement of Return", "Return / Enforcement", Data_test$MaterieAggregated)
Data_test$MaterieAggregated <- gsub("\\&", "and", Data_test$MaterieAggregated)
Data_test$MaterieAggregated <- gsub(" (Asylum)", "", Data_test$MaterieAggregated, fixed = TRUE)
Data_test$MaterieAggregated <- as.factor(Data_test$MaterieAggregated)
Data_test$MaterieAggregated <- relevel(Data_test$MaterieAggregated, ref = "Asylum and Return")

rm(materie_agg_codes)


## Judges anywhere on panel
Data_test$grj1 <- NA # grant rate all other cases panel chair
Data_test$grj2 <- NA # grant rate all other cases second judge
Data_test$grj3 <- NA # grant rate all other cases third judge
Data_test$grp <- NA # grant rate all other cases panel median (grant rate all other cases)

for(i in 1:nrow(Data_test)){ # for each case
  for(b in 1:3){ # and each judge
    a <- Data_test[i, paste0("J", b)] # this is the chair/second/third judge (b) in the ith case
    grj <- data.frame(outcome = Data_test$granted[-i], jop = apply(X=Data_test[-i,c("J1", "J2", "J3")], 1, find_judge)) # outcome in all cases except the ith case and whether the bth judge on the ith case was on the panel or not
    Data_test[i, paste0("grj", b)] <- mean(grj$outcome[grj$jop>0]) # grant rate in all cases the bth judge was on the panel except the ith case
  }
}

Data_test$grp <- apply(X = Data_test[,c("grj1", "grj2", "grj3")], MARGIN = 1, median) # what is the median "leniency" of all the judges on the panels considering all the cases that they were involved in with the exception of the one we're looking at

# regress median panel leniency on case strength indicators, controlling for COO
ols_all_panel <- lm(grp ~ Land + MaterieAggregated + lawyer, data=Data_test)
ols_coo_panel <- lm(grp ~ Land, data=Data_test)

fstat_panel <- anova(ols_all_panel, ols_coo_panel)


## Chair as chair
Data_test$grj2 <- NULL
Data_test$grj3 <- NULL
Data_test$grj1 <- NA

for(i in 1:nrow(Data_test)){
  a <- Data_test[i, "J1"]
  grj <- data.frame(outcome = Data_test$granted[-i], jop = sapply(X=Data_test[-i,"J1"], FUN=find_judge))
  Data_test[i, "grj1"] <- mean(grj$outcome[grj$jop>0])
}

ols_all_chair <- lm(grj1 ~ Land + MaterieAggregated + lawyer, data=Data_test)
ols_coo_chair <- lm(grj1 ~ Land, data=Data_test)

fstat_chair <- anova(ols_all_chair, ols_coo_chair)

# placebo test
placebo_all <- lm(granted ~ Land + MaterieAggregated + lawyer, data=Data_test)
placebo_coo <- lm(granted ~ Land, data=Data_test)

fstat_placebo <- anova(placebo_all, placebo_coo)

table_note <- paste0("Table shows ordinary least squares regressions of the binary appeal outcome (=1 if the appeal is granted) on case characteristics (legal category and legal representation) in Model~1; of the average grant rate of the chair judge across all decisions except case $j$ on case characteristics in Model~2; and of the grant rate of the median judge on the panel across all decisions except case $j$ on case characteristics in Model~3. All models control for country of origin fixed effects. The baseline legal category is `asylum and return'. RR indicates reconsideration requests (following initial rejections).  The joint $F$-test reports the $p$-value from the null hypothesis that all included case strength characteristics are not predictive of the outcome. The sample consists of all ", nrow(Data_test),  " published three-judge panel decisions (granted/rejected) on cases submitted in 2007.")

cov_labels <- gsub("MaterieAggregated", "", attr(coefficients(placebo_all), "names")[-grep("^Land|(Intercept)|Constant", attr(coefficients(placebo_all), "names"))])
cov_labels <- gsub("lawyerTrue", "Lawyer/paralegal", cov_labels)


# table 2

table <- capture.output(stargazer(placebo_all, ols_all_chair, ols_all_panel, 
                                  covariate.labels = cov_labels, 
                                  dep.var.labels = c("Appeal granted", "Leniency Chair", "Leniency Median"),
          add.lines = list(c(paste0("Country FE (#", length(unique(Data_test$Land))-1, ")"), "yes", "yes", "yes"), 
                           c("Joint Test p-value", round3(fstat_placebo$`Pr(>F)`[2]), 
                             round3(fstat_chair$`Pr(>F)`[2]), round3(fstat_panel$`Pr(>F)`[2]))), omit = c(grep("^Land", attr(coefficients(placebo_all), "names"), value=TRUE), "(Intercept)", "Constant"), omit.stat = c("ser","adj.rsq", "f"), 
          notes.append = TRUE, notes=table_note,
          out=paste0("RandomizationTest.tex"),
          font.size="tiny", single.row = TRUE))

rm(Data_test)


############################################
############### ANALYSIS ###################
############################################


###############
# Definitions #
###############

posteriors <- ls(.GlobalEnv, pattern = "posterior_")

model <- gsub("^posterior_", "", posteriors)
model_p <- grep("_p$", model, value = TRUE)
model_indep <- grep("_independent$", model, value = TRUE)
model_s <- grep("_senior$", model, value = TRUE)
model <-
  model[!(model %in% model_p) &
          !(model %in% model_indep) &
          !(model %in% model_s) & model != "mean"]

n_models <- length(model)
n_judges <- length(JudgeName)

## containers
irs <- matrix(NA, nrow = 3, ncol = 3)
lls <- DICs <- dfs <- rep(NA, n_models)
lls_indep <- rep(NA, 3)

# names
names(lls) <-
  names(DICs)  <-
  names(dfs) <- model

rownames(irs) <- names(lls_indep) <- c("mixture", "chair", "median")
colnames(irs) <- c("Inconsistency Rate", "CI_LB95", "CI_UB95")



## extract preference estimates etc from models

for (obj in model) {
  place <- grep(obj, model)
  posterior <-  eval(parse(text = paste0("posterior_", obj)))
  
  ## extract the ll
  lls[place] <- ll.est <- mean(rstan::extract(posterior, "ll")$ll)
  ll.var <- var(rstan::extract(posterior, "ll")$ll)
  DICs[place] <-
    DIC <-
    -2 * ll.est + ll.var / 2
  
  ## extract the number of parameters
  p_mixture <- ifelse(obj == "mixture", 1, 0)
  p_countries <- posterior@par_dims$phi
  p_judges <-
    ifelse(obj != "null", posterior@par_dims$theta, 1) # here's the constant for the null
  dfs[place] <- df <- (p_mixture + p_countries + p_judges)
  
 
  if (obj == "mixture") {
    # extract delta (lambda)
    delta.chain <- rstan::extract(posterior, "delta")$delta
    delta.est <- mean(delta.chain)
    delta.se <- sd(delta.chain)
    delta.ci <- quantile(delta.chain, c(0.025, 0.975))
    print(paste0("Estimated mixture parameter lambda_1 is ", round(delta.est,2), " (95% CI ", round(delta.ci[1],2), "-", round(delta.ci[2],2), ")"))
    
    # get (counterfactual) inconsistency rates
    inconsist.chain <- rstan::extract(posterior, "inconsist")$inconsist
    inconsist.est <- mean(inconsist.chain)
    inconsist.ci <- quantile(inconsist.chain, c(0.025, 0.975))
    
    inconsist.m.chain <- rstan::extract(posterior, "inconsist_median")$inconsist_median
    inconsist.m.est <- mean(inconsist.m.chain)
    inconsist.m.ci <- quantile(inconsist.m.chain, c(0.025, 0.975))
    
    inconsist.c.chain <- rstan::extract(posterior, "inconsist_chair")$inconsist_chair
    inconsist.c.est <- mean(inconsist.c.chain)
    inconsist.c.ci <- quantile(inconsist.c.chain, c(0.025, 0.975))
    
    inconsist.chain.all <-
      cbind(inconsist.chain, inconsist.c.chain, inconsist.m.chain)
    inconsist.est.all <-
      c(inconsist.est, inconsist.c.est, inconsist.m.est)
    inconsist.ci.all <-
      rbind(inconsist.ci, inconsist.c.ci, inconsist.m.ci)
    
    probs.inconsist <-
      matrix(NA,
             ncol = length(inconsist.est.all),
             nrow = length(inconsist.est.all))
    colnames(probs.inconsist) <-
      rownames(probs.inconsist) <-
      names(inconsist.est.all) <-
      rownames(inconsist.ci.all) <-
      colnames(inconsist.chain.all) <-
      c("mixture", "always_chair", "always_median")
    
    
    for (i in 1:length(inconsist.est.all)) {
      for (b in 1:length(inconsist.est.all)) {
        probs.inconsist[b, i] <-
          round(mean(inconsist.chain.all[, b] > inconsist.chain.all[, i]), 3)
      }
    }
    for (i in 1:length(inconsist.est.all)) {
      probs.inconsist[i, i] <- NA
    }
    
    print(xtable(data.frame(
      cbind(inconsist.est.all, inconsist.ci.all)
    ), digits = c(0, 3, 3, 3)),
    file = "CounterfactualIR.tex")
    print(xtable(probs.inconsist), file = "ProbLargerCounterfactualIR.tex")
    
    print(paste0("% increase in inconsistency rate if chair was able to dictate all outcomes: ", round(((inconsist.c.est/inconsist.est)-1)*100,1)))
    
  }
  
  if (obj %in% c("mixture", "chair", "median")) {
   
    ## extract the preference estimates
    theta.chain <- rstan::extract(posterior, "theta")$theta
    
    theta.est <- apply(theta.chain, c(2), mean)
    theta.ci.95 <-
      t(apply(theta.chain, c(2), quantile, c(0.025, 0.975)))
    theta.ci.90 <- t(apply(theta.chain, c(2), quantile, c(0.05, 0.95)))
    names(theta.est) <-
      rownames(theta.ci.95) <- rownames(theta.ci.90) <- JudgeName
    
    probs <- matrix(NA, ncol = n_judges, nrow = n_judges)
    colnames(probs) <-
      rownames(probs) <- colnames(theta.chain) <- JudgeName
    for (i in JudgeName) {
      for (b in JudgeName) {
        probs[b, i] <- round(mean(theta.chain[, b] > theta.chain[, i]), 3)
      }
    }
    for (i in 1:n_judges) {
      probs[i, i] <- NA
    }
    sortorder <- order(theta.est)
    probs <-
      probs[rev(sortorder), rev(sortorder)] # sort by pref estimate
    
    print(xtable(probs),
          file = paste0("ProbJudgeDifferences_", obj, ".tex"))
    
    # extract party-specific pref estimates
    posterior_p <- eval(parse(text = paste0("posterior_", obj, "_p")))
    theta.est.p <-
      apply(rstan::extract(posterior_p, "theta")$theta, c(2), mean)
    theta.ci.95.p <-
      t(apply(
        rstan::extract(posterior_p, "theta")$theta,
        c(2),
        quantile,
        c(0.025, 0.975)
      ))
    theta.ci.90.p <-
      t(apply(
        rstan::extract(posterior_p, "theta")$theta,
        c(2),
        quantile,
        c(0.05, 0.95)
      ))
    theta.chain.p <- rstan::extract(posterior_p, "theta")$theta
    names(theta.est.p) <-
      rownames(theta.ci.95.p) <-
      rownames(theta.ci.90.p) <- colnames(theta.chain.p) <- PartyName
    
    sortorder.p <- order(theta.est.p)
    n_parties <- length(PartyName)
    
    # table whether party-specific pref estimates are different from each other (one-sided tests)
    probs.p <- matrix(NA, ncol = n_parties, nrow = n_parties)
    colnames(probs.p) <- rownames(probs.p) <- PartyName
    for (i in PartyName) {
      for (b in PartyName) {
        probs.p[b, i] <-
          round(mean(theta.chain.p[, b] > theta.chain.p[, i]), 3)
      }
    }
    for (i in 1:n_parties) {
      probs.p[i, i] <- NA
    }
    probs.p <- probs.p[rev(sortorder.p), rev(sortorder.p)]
    print(xtable(probs.p),
          file = paste0("ProbMoreLenientByParty_", obj, ".tex"))
    
    ## create plots for three best fitting models with judge- and party-specific pref estimates side by side

    JudgeColor <- apply(JData[, c("red", "green", "blue")], 1,
                        function(x)
                          rgb(
                            red = x[1],
                            green = x[2],
                            blue = x[3],
                            maxColorValue = 255
                          ))
    
    plotname <- paste0("JudgePreferences_", obj, ".png")
    png(file = plotname, width=12, height=6, units="in", res=600)
    
    layout(
      mat = matrix(c(1, 2), 1),
      width = c(6, 6),
      height = c(7, 7)
    )
    par(mar = c(4, 4, 2, 0.5))
    
    y_axis <- "Preferred Fraction of Appeals Granted"
    
    plot(
      theta.est[sortorder],
      1:n_judges,
      xlim = c(0, 0.66),
      axes = FALSE,
      pch = 16,
      xlab = "Judge Preferences (Posterior Means)",
      ylab = "",
      col = JudgeColor[sortorder],
      cex = 1.5,
      cex.lab = 0.8
    )
    
    for (i in 1:n_judges)
      abline(h = i, col = "light grey", lty = 3)
    for (i in 1:n_judges)
      lines(
        c(theta.ci.95[sortorder[i], 1], theta.ci.95[sortorder[i], 2]),
        rep(i, 2),
        lwd = 1,
        col = (JudgeColor[sortorder])[i]
      )
    for (i in 1:n_judges)
      lines(
        c(theta.ci.90[sortorder[i], 1], theta.ci.90[sortorder[i], 2]),
        rep(i, 2),
        lwd = 3,
        col = (JudgeColor[sortorder])[i]
      )
    abline(v = mean(theta.est), lty = 2)
    axis(
      2,
      at = 1:n_judges,
      labels = names(theta.est)[sortorder],
      las = 1,
      cex.axis = 0.75
    )
    
    PartyColor <- unique(JudgeColor)
    legend(
      "bottomright",
      legend = PartyName[sortorder.p],
      col = PartyColor[sortorder.p],
      lwd = 2,
      pch = 16,
      cex = 0.8,
      bg = "white"
    )
    ppts <- seq(0, 1, 0.1)
    axis(1,
         at = ppts,
         labels = ppts,
         cex.axis = 0.8)
    axis(3,
         at = ppts,
         labels = ppts,
         cex.axis = 0.8)
    
    # party plot
    y_upper <- max(theta.ci.95.p + 0.03)
    x_parties <- c(NA, theta.est.p[sortorder.p], NA)
    y_parties <- 0:(n_parties + 1)
    
    par(mar = c(4, 0.5, 2, 4))
    
    plot(
      x_parties,
      y_parties,
      xlim = c(0, y_upper),
      axes = FALSE,
      pch = 16,
      xlab = "Judge Preferences by Party (Posterior Means)",
      ylab = "",
      col = c(NA, PartyColor[sortorder.p], NA),
      cex = 1.5,
      cex.lab = 0.8
    )
    
    for (i in 1:n_parties) {
      abline(h = i, col = "light grey", lty = 3)
      lines(
        c(theta.ci.95.p[sortorder.p[i], 1], theta.ci.95.p[sortorder.p[i], 2]), 
        rep(i, 2),
        lwd = 1,
        col = (PartyColor[sortorder.p])[i]
      )
      lines(
        c(theta.ci.90.p[sortorder.p[i], 1], theta.ci.90.p[sortorder.p[i], 2]),
        rep(i, 2),
        lwd = 3,
        col = (PartyColor[sortorder.p])[i]
      )
    }
    
    abline(v = mean(theta.est), lty = 2)
    
    
    # neighbouring parties different from each other?
    ps <- PartyName[sortorder.p]
    
    # sp vs cvp
    segments(
      x0 = theta.ci.95.p[ps[n_parties], 2] + 0.01,
      x1 = theta.ci.95.p[ps[n_parties], 2] + 0.02,
      y1 = which(ps == ps[n_parties]),
      y0 = which(ps == ps[n_parties])
    )
    segments(
      x0 = theta.ci.95.p[ps[n_parties], 2] + 0.01,
      x1 = theta.ci.95.p[ps[n_parties], 2] + 0.02,
      y1 = which(ps == ps[n_parties - 1]),
      y0 = which(ps == ps[n_parties - 1])
    )
    segments(
      x0 = theta.ci.95.p[ps[n_parties], 2] + 0.02,
      x1 = theta.ci.95.p[ps[n_parties], 2] + 0.02,
      y1 = which(ps == ps[n_parties - 1]),
      y0 = which(ps == ps[n_parties])
    )
    segments(
      x0 = theta.ci.95.p[ps[n_parties], 2] + 0.02,
      x1 = theta.ci.95.p[ps[n_parties], 2] + 0.015,
      y1 = mean(c(which(ps == ps[n_parties - 1]), which(ps == ps[n_parties]))),
      y0 = mean(c(which(ps == ps[n_parties - 1]), which(ps == ps[n_parties])))
    )
    
    text(
      x = theta.ci.95.p[ps[n_parties], 2] + 0.01,
      y = mean(c(which(ps == ps[n_parties - 1]), which(ps == ps[n_parties]))),
      label = paste0("Pr(", ps[n_parties], " > ", ps[n_parties - 1], ") = ", 
                     round(probs.p[ps[n_parties], ps[n_parties - 1]], 2)),
      cex = 0.75,
      adj = 1
    )
    
    # cvp indep line:
    segments(
      x0 = theta.ci.95.p[ps[n_parties - 1], 2] + 0.01,
      x1 = theta.ci.95.p[ps[n_parties - 1], 2] + 0.02,
      y1 = which(ps == ps[n_parties - 2]),
      y0 = which(ps == ps[n_parties - 2])
    )
    segments(
      x0 = theta.ci.95.p[ps[n_parties - 1], 2] + 0.01,
      x1 = theta.ci.95.p[ps[n_parties - 1], 2] + 0.02,
      y1 = which(ps == ps[n_parties - 1]),
      y0 = which(ps == ps[n_parties - 1])
    )
    segments(
      x0 = theta.ci.95.p[ps[n_parties - 1], 2] + 0.02,
      x1 = theta.ci.95.p[ps[n_parties - 1], 2] + 0.02,
      y1 = which(ps == ps[n_parties - 1]),
      y0 = which(ps == ps[n_parties - 2])
    )
    segments(
      x0 = theta.ci.95.p[ps[n_parties - 1], 2] + 0.02,
      x1 = theta.ci.95.p[ps[n_parties - 1], 2] + 0.015,
      y1 = mean(c(which(ps == ps[n_parties - 1]), which(ps == ps[n_parties - 2]))),
      y0 = mean(c(which(ps == ps[n_parties - 1]), which(ps == ps[n_parties - 2])))
    )
    
    text(
      x = theta.ci.95.p[ps[n_parties - 1], 2] + 0.01,
      y = mean(c(which(ps == ps[n_parties - 1]), which(ps == ps[n_parties -  2]))),
      label = paste0("Pr(", ps[n_parties - 1], " > ", ps[n_parties - 2], ") = ", 
                     round(probs.p[ps[n_parties - 1], ps[n_parties - 2]], 2)),
      cex = 0.75,
      adj = 1
    )
    
    # cvp fdp line:
    segments(
      x0 = theta.ci.95.p[ps[n_parties - 2], 2] + 0.01,
      x1 = theta.ci.95.p[ps[n_parties - 2], 2] + 0.02,
      y1 = which(ps == ps[n_parties - 2]),
      y0 = which(ps == ps[n_parties - 2])
    )
    segments(
      x0 = theta.ci.95.p[ps[n_parties - 2], 2] + 0.01,
      x1 = theta.ci.95.p[ps[n_parties - 2], 2] + 0.02,
      y1 = which(ps == ps[n_parties - 3]),
      y0 = which(ps == ps[n_parties - 3])
    )
    segments(
      x0 = theta.ci.95.p[ps[n_parties - 2], 2] + 0.02,
      x1 = theta.ci.95.p[ps[n_parties - 2], 2] + 0.02,
      y1 = which(ps == ps[n_parties - 3]),
      y0 = which(ps == ps[n_parties - 2])
    )
    segments(
      x0 = theta.ci.95.p[ps[n_parties - 2], 2] + 0.02,
      x1 = theta.ci.95.p[ps[n_parties - 2], 2] + 0.025,
      y1 = mean(c(which(ps == ps[n_parties - 3]), which(ps == ps[n_parties - 2]))),
      y0 = mean(c(which(ps == ps[n_parties - 3]), which(ps == ps[n_parties - 2])))
    )
    
    text(
      x = theta.ci.95.p[ps[n_parties - 2], 2] + 0.03,
      y = mean(c(which(ps == ps[n_parties - 3]), which(ps == ps[n_parties - 2]))),
      label = paste0("Pr(", ps[n_parties - 2], " > ", 
                     ps[n_parties -  3], ") = ", 
                     round(probs.p[ps[n_parties - 2], ps[n_parties - 3]], 2)),
      cex = 0.75,
      adj = 0
    )
    
    # FDP SVP line:
    segments(
      x0 = theta.ci.95.p[ps[n_parties - 3], 2] + 0.01,
      x1 = theta.ci.95.p[ps[n_parties - 3], 2] + 0.02,
      y1 = which(ps == ps[n_parties - 4]),
      y0 = which(ps == ps[n_parties - 4])
    )
    segments(
      x0 = theta.ci.95.p[ps[n_parties - 3], 2] + 0.01,
      x1 = theta.ci.95.p[ps[n_parties - 3], 2] + 0.02,
      y1 = which(ps == ps[n_parties - 3]),
      y0 = which(ps == ps[n_parties - 3])
    )
    segments(
      x0 = theta.ci.95.p[ps[n_parties - 3], 2] + 0.02,
      x1 = theta.ci.95.p[ps[n_parties - 3], 2] + 0.02,
      y1 = which(ps == ps[n_parties - 3]),
      y0 = which(ps == ps[n_parties - 4])
    )
    segments(
      x0 = theta.ci.95.p[ps[n_parties - 3], 2] + 0.02,
      x1 = theta.ci.95.p[ps[n_parties - 3], 2] + 0.025,
      y1 = mean(c(which(ps == ps[n_parties - 3]), which(ps == ps[n_parties - 4]))),
      y0 = mean(c(which(ps == ps[n_parties - 3]), which(ps == ps[n_parties - 4])))
    )
    text(
      x = theta.ci.95.p[ps[n_parties - 3], 2] + 0.03,
      y = mean(c(which(ps == ps[n_parties - 3]), which(ps == ps[n_parties - 4]))),
      label = paste0("Pr(", ps[n_parties - 3], " > ", ps[n_parties - 4], ") = ", 
                     round(probs.p[ps[n_parties - 3],  ps[n_parties - 4]], 1)),
      cex = 0.75,
      adj = 0
    )
    axis(
      4,
      at = 1:n_parties,
      labels = names(theta.est.p)[sortorder.p],
      las = 1,
      cex.axis = 0.75
    )
    
    ppts <- seq(0, 1, 0.1)
    axis(1,
         at = ppts,
         labels = ppts,
         cex.axis = 0.8)
    axis(3,
         at = ppts,
         labels = ppts,
         cex.axis = 0.8)

    dev.off()
    
    ## Print out judge-specific Preference Estimates for Chair, Median, Mixture
    PreferenceEstimatesModel <-
      data.frame(
        PointEstimate = theta.est[rev(sortorder)],
        CI95_LB = theta.ci.95[rev(sortorder), 1],
        CI95_UB = theta.ci.95[rev(sortorder), 2]
      )
    
    if (obj == "mixture") {
      delta <- c(delta.est, delta.ci)
      names(delta) <- names(PreferenceEstimatesModel)
      PreferenceEstimatesModel <-
        rbind(delta, PreferenceEstimatesModel)
    }
    
    ## Print out party-specific Preference Estimates for Chair, Median, Mixture
    PreferenceEstimatesModelP <-
      data.frame(
        PointEstimate = theta.est.p[rev(sortorder.p)],
        CI95_LB = theta.ci.95.p[rev(sortorder.p), 1],
        CI95_UB = theta.ci.95.p[rev(sortorder.p), 2],
        CI90_LB = theta.ci.90.p[rev(sortorder.p), 1],
        CI90_UB = theta.ci.90.p[rev(sortorder.p), 2]
      )
    
    namePE <-
      paste0("PreferenceEstimates_", obj)
    print(xtable(PreferenceEstimatesModel, digits = c(0, 3, 3, 3)),
          file = paste0(namePE, ".tex"))
    print(xtable(PreferenceEstimatesModelP, digits = c(0, 3, 3, 3, 3, 3)),
          file = paste0(namePE, "_p.tex"))
    
    # average judge preference
    print(paste0("Average Judge Preference (", obj, " model): ", round(mean(theta.est),3)))
    
    
    ## extract the inconsistency rate (based on case-space assumptions)
    irs[obj, 1] <-
      inconsist.est <-
      mean(rstan::extract(posterior, "inconsist")$inconsist)
    irs[obj, 2:3] <-
      inconsist.ci <-
      quantile(rstan::extract(posterior, "inconsist")$inconsist,
               c(0.025, 0.975))
    
    if (obj == "mixture") {
     
      # Fitted Value Distribution Graph in appendix
      pi.lk <- matrix(NA,
                      nrow = dim(extract(posterior, "pi")$pi)[1],
                      ncol = dim(extract(posterior, "pi")$pi)[2])
      colnames(theta.chain) <- JudgeName
      panels <- Data[, c("J1", "J2", "J3")]
      temp <- panels
      
      for (th in 1:nrow(theta.chain)) {
        for (i in colnames(panels)) {
          temp[, i] <-
            theta.chain[th, match(panels[, i], names(theta.chain[th, ]))]
        }
        
        # using reference country (Sri Lanka), not actual COO
        chair <- temp[, 1]
        median <- t(apply(temp, MARGIN = 1, sort))[, 2]
        pi.lk[th, ] <-
          chair * delta.chain[th] + median * (1 - delta.chain[th])
      }
      
      plotname <-
        paste0("FittedValueDistribution_LK_", obj , ".png")
      
      png(plotname, width=4, height=4, units="in", res=600)
      par(mar = c(4, 3, 0, 0) + 0.1, mfrow = c(1, 1))
      
      hist(
        pi.lk,
        br = seq(0, max(pi.lk) + 0.01, 0.01),
        main = "",
        xlab = "",
        ylab = "",
        xlim = c(0, max(pi.lk) + 0.01),
        col = "grey",
        border = NA,
        freq = FALSE,
        cex.lab = 0.7,
        cex.sub = 0.7,
        cex.axis = 0.7
      )
      
      mtext(
        side = 1,
        text = "Fitted Probability",
        line = 2,
        cex = 0.7
      )
      mtext(
        side = 1,
        text = paste("Inconsistency Rate: ", round(inconsist.est, 3), sep = ""),
        line = 3,
        cex = 0.7
      )
      mtext(
        side = 2,
        text = "Density",
        line = 2,
        cex = 0.7
      )
      
      abline(v = mean(pi.lk), lty = 2)
      
      dev.off()
      print(paste0("Mean predicted probability of appeal success (mixture model) holding country of origin constant (at Sri Lanka): ", round(mean(pi.lk),2), "; Mean grant rate: ", round(mean(Data$granted), 2)))
    }
  }
}

# table 3
fittable <-
  data.frame(round(DICs, 1), round(lls, 1), dfs) 
colnames(fittable) <-
  c("DIC", "LL", "Parameters") 
fittable <- fittable[order(fittable$DIC), ]
print(xtable(fittable, digits = c(0, 1, 1, 0)), file = "FitTablePosterior.tex")

# compare lls with lls_indept
lls_indep["mixture"] <-
  mean(rstan::extract(posterior_mixture_independent, "ll")$ll)
lls_indep["chair"] <-
  mean(rstan::extract(posterior_chair_independent, "ll")$ll)
lls_indep["median"] <-
  mean(rstan::extract(posterior_median_independent, "ll")$ll)
lls_table <-
  data.frame(LL = fittable$LL[rownames(fittable) %in% c("mixture", "chair", "median")],
             LL_indep = lls_indep)

print(xtable(lls_table, digits = c(0,1, 1)), file = "LL_comp.tex")

# mean model ll:
print(paste0("mean model ll: ", round(mean(
  rstan::extract(posterior_mean, "ll")$ll
), 1)))

# save inconsistency rate table
print(xtable(irs, digits = c(0,3, 3, 3)), file = "InconsistencyRate.tex")


# inconsistency rate independence assumption
irs_indep.chain <- rstan::extract(posterior_mixture_independent, "inconsist")$inconsist
irs_indep <- round(quantile(irs_indep.chain,c(0.5, 0.025, 0.975))*100,1)
print(paste0("IR mixture model (independence assumption): ", irs_indep[1], "% (95% CI: ", irs_indep[2], "%-", irs_indep[3], "%)"))

### split sample model

posterior <- posterior_mixture_senior

delta1.chain <- rstan::extract(posterior, "delta1")$delta1
delta1.est <- mean(delta1.chain)
delta1.se <- sd(delta1.chain)
delta1.ci <- quantile(delta1.chain, c(0.025, 0.975))
delta0.chain <- rstan::extract(posterior, "delta0")$delta0
delta0.est <- mean(delta0.chain)
delta0.se <- sd(delta0.chain)
delta0.ci <- quantile(delta0.chain, c(0.025, 0.975))

prob.lambda <- mean(delta1.chain > delta0.chain)

print(paste0("mixture parameter for the chair's power (senior chairs): ", round(delta1.est,2), " (95% CI: ", round(delta1.ci[1],2), "-", round(delta1.ci[2],2), ")"))
print(paste0("mixture parameter for the chair's power (junior chairs): ", round(delta0.est,2), " (95% CI: ", round(delta0.ci[1],2), "-", round(delta0.ci[2],2), ")"))

# pref estimates
theta.est <-
  apply(rstan::extract(posterior, "theta")$theta, c(2), mean)
theta.ci.95 <-
  t(apply(
    rstan::extract(posterior, "theta")$theta,
    c(2),
    quantile,
    c(0.025, 0.975)
  ))
theta.ci.90 <-
  t(apply(
    rstan::extract(posterior, "theta")$theta,
    c(2),
    quantile,
    c(0.05, 0.95)
  ))

sortorder <- order(theta.est, decreasing = TRUE)
names(theta.est) <- JudgeName

## Print out Preference Estimates
PreferenceEstimatesModel <-
  data.frame(PointEstimate = theta.est[sortorder],
             CI95_LB = theta.ci.95[sortorder, 1],
             CI95_UB = theta.ci.95[sortorder, 2])
delta1 <- c(delta1.est, delta1.ci)
delta0 <- c(delta0.est, delta0.ci)
delta <- rbind(delta1, delta0)
colnames(delta) <- colnames(PreferenceEstimatesModel)
rownames(delta) <-
  c("Delta Senior Chair", "Delta Junior Chair")
PreferenceEstimatesModel <-
  rbind(delta, PreferenceEstimatesModel)


namePE <-
  "PreferenceEstimatesSenioritySplit_mixture.tex"
print(xtable(PreferenceEstimatesModel, digits = c(0, 3, 3, 3)), file =
        namePE)

print(paste0("Pr(lambda_1 Senior Chair Judge > lambda_1 Junior Chair Judge): ", round3(prob.lambda)))
